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Abstract 

The stable and efficient numerical simulation of viscoelastic flows has been a constant struggle due to the 
High Weissenberg Number Problem. While the stability for macroscopic descriptions could be greatly 
enhanced by the log-conformation method as proposed by Fattal and Kupferman, the application of the 
efficient Newton-Raphson algorithm to the full monolithic system of governing equations, consisting of the 
log-conformation equations and the Navier-Stokes equations, has always posed a problem. In particular, it 
is the formulation of the constitutive equations by means of the spectral decomposition that hinders the 
application of further analytical tools. Therefore, up to now, a fully monolithic approach could only be 
achieved in two dimensions, as, e.g., recently shown in [P. Knechtges, M. Behr, S. Elgeti, Fully-implicit 
log-conformation formulation of constitutive laws, J. Non-Newtonian Fluid Mech. 214 (2014) 78-87]. 

The aim of this paper is to find a generalization of the previously made considerations to three 
dimensions, such that a monolithic Newton-Raphson solver based on the log-conformation formulation 
can be implemented also in this case. The underlying idea is analogous to the two-dimensional case, to 
replace the eigenvalue decomposition in the constitutive equation by an analytically more "well-behaved" 
term and to rely on the eigenvalue decomposition only for the actual computation. Furthermore, in order to 
demonstrate the practicality of the proposed method, numerical results of the newly derived formulation are 
presented in the case of the sedimenting sphere and ellipsoid benchmarks for the Oldroyd-B and Giesekus 
models. It is found that the expected quadratic convergence of Newton's method can be achieved. 

Keywords: Log-conformation, Oldroyd-B model, Giesekus model. Finite element method 
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1. Introduction 

Viscoelastic flows are ubiquitous in modem industrial applications. They are essential for the correct 
description of the flow properties of blood, as well as polymer melts, which makes a good understanding of 
the used models necessary for applications ranging from the design of blood pumps [1] to the layout of 
extrusion dies in plastics manufacturing [2]. 

Considering the demands stemming from the non-linear behavior of most of the used models and, at the 
same time, the possibilities given through the advent of the computer age, it has become more and more 
common not to base the model analysis solely on pure analytic grounds, but also to perform numerical 
simulations, which can be applied to almost arbitrary geometries and domains. In the past, the macroscopic 
descriptions have been quite dominant, whereas micro-macro simulations based on stochastic differential 
equations are now gaining importance [3, 4]. Although the latter offer a greater flexibility with respect 
to the modeling of the underlying molecular dynamics, the former are still quite popular due to their 
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lower computational cost. Since this is important for the application of fhe numerical mefhods fo complex 
geometries, this paper seeks a description in the macroscopic framework. More specifically, we will consider 
fhe Oldroyd-B [5] and fhe Giesekus model [6]. The applicabilify of our mefhods, however, is nof limifed fo 
fhese two models. 

Simulfaneously fo fhe advenf of numerical mefhods in fhe simulafion of viscoelasfic models, fhe High 
Weissenberg Number Problem arose [7, 3]. Wifh fhe Weissenberg number being a dimensionless consfanf 
fhaf weighfs fhe confribufion of fhe viscoelasfic equafions fo fhe descripfion of fhe full sysfem, fhis absfracf 
ferm expresses fhe empirical facf fhaf, wifh increasing Weissenberg number, numerical simulations fend fo 
fail. In facf, fhe range of affainable Weissenberg numbers fumed ouf fo be quife offen lower fhan whaf was 
measured in experimenfs, fhus reducing fhe predicfive power of simulations. 

The most recent and quite successful approaches tackling the High Weissenberg Number Problem are the 
log-conformation methods, first considered in [8] in order to better resolve exponential stress-boundary layers 
in regions of high sfrain. Alfhough fhey do nof solve fhe High Weissenberg Number Problem complefely, 
they address the subproblem that numerical simulations do not necessarily preserve the positive-definiteness 
of fhe conformafion fensor; a properfy always fulfilled by fhe undiscrefized equafions [9]. The laffer was 
found fo be crucial for a numerical simulafion nof fo fail. The underlying idea of fhe log-conformafion 
mefhods is as simple as if is powerful: The so far primal degree of freedom — fhe conformafion fensor cr — 
is replaced by its logarithm T. Hence, cr is obtained by means of fhe mafrix-exponenfial funcfion exp 'P, 
which aufomafically ensures fhaf cr remains posifive-definife. 

This, however, comes af fhe cost of finding a suifable replacemenf for fhe corresponding consfifufive 
equation. The way fhe original mefhod [8] pursues is rafher unusual, compared fo ofher partial differential 
equafions, in the requirement of an eigenvalue decomposition of T. In particular, if is this spectral 
decomposition that hinders the direct application of fhe Newton-Raphson algorifhm to the full set of 
partial differential equations. More specifically, fhe Newfon-Raphson mefhod involves a determination 
of derivafives wifh respect to the 'P degrees of freedom, including the derivatives of fhe eigenvalues and 
eigenvectors. Nevertheless, considering the derivatives of eigenvectors, it is known that they become 
singular in the case of degenerate eigenvalues due to the ambiguity in the eigenvectors. As a remedy for this 
and for fhe difficulfy of faking fhe derivative of fhe mafrix-exponenfial funcfion, firsf affempfs resorfed fo 
fhe approximafion of fhe Jacobian mafrix by difference quofienfs [10,11]. 

Even though a first analysis was conducted for the two-dimensional Leonov model in [12], it was not until 
the work in [13] and [14] that the Jacobian matrix was derived by pure analytic means in two dimensions for 
a broader class of models. As a continuation of these earlier works, this paper is devoted to a generalization 
to arbitrary dimensions, along which we will also bridge the gap between the two expositions in [13] and 
[14]. 

In particular, we will not only discuss the derivation of a new consfifufive equation in fhe firsf section, 
but also describe the numerical implications in the case of an implementation into an existing Galerkin/Least- 
Squares (GLS) Navier-Stokes solver in the succeeding section. The results of fhis solver are subsequenfly 
used in Secfion 4 fo sfudy fhe falling sphere benchmark, where a sphere of radius R sediments along the 
centerline of a fube of radius 2R. In order fo demonsfrafe fhe applicabilify fo fruly fhree-dimensional flows, 
a modification of fhe same benchmark wifh a fri-axial ellipsoid is considered as well. 

Alfhough our mofivafion stems mostly from the numerical side, the proposed equations are purely 
analytic and as such may also serve as a new tool in future analytic studies; to the author's knowledge, this 
is the first time that the constitutive equations for the log-conformation formulation can be stated in a closed 
form in fhis generalify. 

2. Theory 

The aim of fhis secfion is fhe derivation of an alfemafive consfifufive equation wifh T as a new primal 
variable. Sfarfing point is the original constitutive equation in terms of fhe conformafion fensor cr and 
fhe velocity u. Both are fields fhat, given boundary and inifial condifions, have to be determined over a 
time-span [0, T] and a d-dimensional domain Q c IR*^. 
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Following the exposition in [13], we consider constitutive models of the form 

1 

dtcr + (m • V)cr + [cr, Q(m)] - e{u)cr — cre{u) = —rf’(£r), (1) 

A 

where e{u) = 1 ^Vm + denotes the strain tensor, Q(m) = 1 ^Vm - the vorticity tensor, A the relaxation 
time, and P{o-) an analytic function. The bracket [cr, Q(m)] is the so-called commutator, which is defined as 

[cr, Q(m)] = (tQ.{u) - Q(H)cr. 

Common choices for P(cr) are f’(cr) = cr - 1, leading fo fhe Oldroyd-B model [5], or P(cr) = cr - 1 -t- a(cr - 1)^ 
with a e [0,1] in the Giesekus model [6]. Generalizations of the subsuming methods to the Johnson-Segalman 
model, as, e.g., done in [14], or other models are in principle possible, but omitted here for the sake of brevity. 

Since the velocity field u is nof defermined so far, we have fo combine fhe consfifufive equations wifh 
fhe Navier-Stokes equafions in order fo obfain a complete system of parfial differenfial equations. More 
specifically, fhe Navier-Sfokes equafions are given by 


V • M = 0 

Up (2) 

p(dt-hM-V)M-hVp-2psV-e(M)- yV-(o--l) = 0, 

wifh densify p, as well as solvenf and polymeric viscosify consfanfs ps and pp, respectively. 

Furfhermore, we will infroduce function spaces 'H and 'H', which for fhe momenf could be chosen as 

perfecfly smoofh, i.e., “K = “K' = C“([0, T] X Q), and fhe derived spaces 

H = Hsym = {XeH\X'^ = X], 

H' = = [Z e H' I = Z]. 

The cenfral sfafemenf of fhis paper fhen reads: 

Theorem 1. Let the velocity field u be given with e{u) e and 0(m) e H'. If *P e Hsym satisfies 

+ (« . V)^ + [‘P,Q(«)] - ^ J J/(z - z')^e{u)^dzdz' = -\p(e'^) e-''’ (3) 


with 


2x X 

fix) ^ _ tanh(x/2)' 

and F chosen as a closed path surrounding the spectrum o/T in {z eC\\ Im(z)| < ti], then cr = exp *F e Hsy„ solves 
the original constitutive equation (1). 

Before we come to the proof of this theorem, we need to consider certain properties of the relevant terms. 
The first thing to notice in Eq. (3) is the double integral, which is similar to the familiar Cauchy infegral 
from complex analysis. One of fhe main differences fo fhe ordinary Cauchy infegral, however, is fhaf fhe 
scalar ratio has been replaced by the resolvent l/(z - 'F) (zl - where 1 is the identity matrix and 
indicates the matrix inverse. For smooth function spaces, it can be deduced that, at a specific insfanf of 
space and time, the resolvent exists if and only if z does not equal any of the eigenvalues of *P(f, x), which are 
all real-valued. Encircling these poles with our integration path F subsequently gives us, by the same means 
as in the Cauchy integral setting, some information on / at these poles, but with the additional complexity 
that we have to deal with the matrix algebra. 

Although we will further use this idea of numerically evaluating fhe infegrals af fhe eigenvalues lafer 
on, we will now leave fhe seffing of smoofh function spaces. Insfead, we consider fhe more general case of 
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Re(z) 


Figure 1: Illustration of a particular choice of the integration path F, as used in Theorem 1. Here, choosing F as an ellipse, with 
semi-major axis greater than ||'P||h and semi-minor axis smaller than n, ensures that the spectrum of T is enclosed by F, while the poles 
of /, especially ±2ni, do not contribute to the integral. 


choosing 9^' as a Banach space and 99 c 99' as a commutative Banach algebra. This opens up the door to a 
variety of spaces as they are used in the analysis of parfial differenfial equations. An example sef of spaces 
in fhis more general selling would be fhe Sobolev-based spaces 

99 = Ci([0,T],H*-i(Q)) n C°([0,T],H*(O)), 

99' = C°([0,T],H®-i(Q)), 


with s > dfl and Q c ]R‘^ being a Lipschitz-bounded domain [15]. It should be stressed that the mathematical 
discussion here is not limited to these spaces, and for the general requirements on 99 and 99' we refer to the 
appendix of [13]. 

Considering whefher on these general spaces the double integral is well defined, fhe set of values for 
which fhe resolvenf is nof defined is no longer resfricted fo fhe disfrncf eigenvalues, but may in fact be 
larger, although still real-valued. Therefore, in order fo separafe fhe ferminology from fhe mafrix algebra, 
fhis sef is called specfrum in fhe general Banach algebra setting. The generalizafion of Cauchy's infegral fo 
fhe fheory of Banach algebras^ is known as Dunford's infegral and if is fhe essenfial ingredienf fo define a 
functional calculus of holomorphic functions on fhese algebras [16,17]. More specifically, for a function g 
fhaf is holomorphic in fhe neighborhood of fhe specfrum of *F, fhe H-valued funcfion is defined as 




1 r g(z) 

2ni Jy z 


dz, 


(4) 


where T is a contour surrounding the spectrum of T wifhin the same neighborhood. Here, as well as 
in Theorem 1, it is assumed that T encircles the spectrum only once. An immediate consequence of this 
definition is that it allows us to explain the exponential function of T, which also has to be an element of 
hfsym* 

A discrepancy between the integral in Theorem 1 and the usual Dunford integral is that the argument of 
/ depends on the difference between the two integral variables z and z'. The latter is also what makes it more 
difficult to ensure that the poles of /, especially ±2ni, do not contribute to the integral. In the formulation of 
Theorem 1 this has been realized by restricting the imaginary part of fhe integration path T to the region 
I Imz| < n. An example of a closed curve T fulfilling fhe aforemenfioned criferia is depicfed in Fig. 1, where 


'Throughout this paper we deliberately use the same symbol for the Banach algebra as for its complexification. 
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the fact is also used that the spectrum of *P is always confained in fhe inferval [-||*F||h/ 11‘1'IIh]- 

Despite this minor restriction, most of the properties of Dunford's infegral carry over to the double 
integral as well. One of the more important features is the independence of the integral on the exact contour 
of r, which is a key consequence of Cauchy's fheorem [16, Theorem 3.31]. 

The final reason for not including the poles of / in the integral in Theorem 1 is that we want to express / 
by a Taylor series in the course of fhe proof. 

Lemma 1. Let g be a holomorphic function on a convex domain Q' c C. Moreover, let the ball of radius r, 6,(0), 
be contained in O'. The Taylor series on this ball shall be given by g{z) = L^o Then for every A e H with 
||A||h < r/2, B e H', and T c ^O' a contour around the spectrum of A, it holds 


oo 

= J^b„{A,B}„, 


-dzdz' 


where {A,B}„ denotes the n-th iterated commutator 


{A,B]n -[A,{A,B}„_i] = 2]^ r (-l)'A”-'BAh 

Proof. Wifhouf loss of generalify, Cauchy's fheorem allows us fo choose a confour T wifhin B,/ 2 ( 0 ) fhaf sfill 
surrounds the spectrum of A. Since the Taylor series converges uniformly on every compact subset of 
6,(0), and especially on T - T c 6,(0), one can deduce by similar means as for Dunford's infegral (cf. [16, 
Theorem 10.27]) fhaf g in F(A, B) can be approximafed by fhe Taylor series fo yield an arbifrarily accurafe 
approximation of F{A,B). As such, we can assume g{z - z') = (z - z')". Furfhermore, using a binomial 
expansion (z - z')” = Xl"=o (")(-l)'z"“‘z'* we obfain 




-dzdz' 




-dz' , 


which fogefher wifh (4), or more rigorously [16, Lemma 10.24], yields fhe desired iferafed commufafor 

= ^(”)(-1)'A”-'BA* = {A,B)„. 


1=0 


□ 

This lemma is already sufficienf fo prove Theorem 1 in the case ||*P||h < n, as can be seen by choosing 
g{z) - f{z) = 2 Yj7=o well as Q' = IR + i{-2n, 2n) and r - 2n. Then it becomes apparent that if *P 

fulfills Eq. (3), if also has fo fulfill 

+ {u ■ V)‘P + [T, 0(h)] - 2 2_, e{u)] 2 „ = -jP (e'^) e”"', 

n=0 ' '■ 

where 62 „ denote the even Bernoulli numbers. This is exactly the equation for which the conclusion in 
Theorem 1 has been proven in [13, Theorem 1]. 

The generalization to H'PHh > n is part of fhe 
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Proof of Theorem 1. Assuming T solves Eq. (3), in a first step the Wilcox Lemma [18] is applied to handle the 
derivatives <9f + m • V of the exponential mapping cr = exp T, such that subsequently (3) can be inserted: 


(<9f + M • V) cr 


= f ((<9t + H • V) T) da 

Jo 

= -l j^'e(i-«)'i'p(e'r)e-V'^da 


Jo 


^ u 

f 


The integral involving P is the easiest to handle since all involved terms commute with each other, 

resulting in the contribution -jP{cr). The vorticity term can be simplified by the fundamental theorem of 
calculus, which yields 


- J" e^^ □(«)]£"'*’da = J" da = -e'^D(u) + D(u)e'^ = -[cr, Q(m)] . 

Hence, for cr to fulfill the original constitutive equation it is left to prove 


Rather than directly proving this equality we will follow the argumentation of Theorem 2 in [13], and 
consider a slightly more general equation where *F is replaced by f'P, such that an analytic continuation 
argument in f can be used to bridge the gap to the case ||'P||h < ri. In particular, without loss of generality 
we will assume that T, in addition to the spectrum of 'V, also encloses B7i/2(0), as, e.g., depicted in Fig. 1. Our 
assertion now reads that 


shall hold with 


f £(H))e“^'^ da =e(u)ef^'^ + eP'^e(u) 

Jo 


(5) 


for every jS in a sufficiently small simply-connected neighborhood D of [0,1] U B7 i/(2||>p||h)(0) ^ C. 
It is clear that the right-hand side of Eq. (5) is holomorphic for all jS, due to 




Additionally, using 


111 


z-f'P z-f'V z-jS*P' 
it is evident that the left-hand side is holomorphic for jS e D. 


( 6 ) 
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Restricting ourselves for a moment to |jS| < 7 t/(2||*P||h), we deduce from Lemma 1 fhaf 


J r>i ( °° 

0 -'o [ ^ 


P2n 

(2m)! 


{jS'P,e(H)} 2 „ 




which is essentially the form for which the identity (5) has already been proven in the proof of Theorem 2 in 
[13]. Thus, we are leff wifh applying fhe monodromy fheorem fhat asserts the uniqueness of the analytic 
continuation on D. Thereby, (5) has in particular to hold for jS = 1, and cr solves fhe original consfifufive 
equafion ( 1 ). □ 


3. Numerical implementation 

Given the newly derived constitutive equation (3), we are going to discuss the numerical discretization 
in conjunction with the Finite Element Method (FEM). The first part of fhis secfion will be cenfered around 
fhe formulation of fhe discrefized weak form in ferms of space-fime elemenfs. The second parf will fhen be 
concerned wifh fhe linearizafion of fhe discrefized weak form by means of fhe Newfon-Raphson mefhod. 
In particular, it will also deal with the evaluation of the double integral and its derivatives. It should be 
noted that the two subsections are only loosely coupled and that the discussion of the latter subsection is not 
limited to the discretization scheme we have chosen, but may be easily generalized to other schemes. 

3.1. Discretization 

As in the preceding paper [13], we will use a mixture of a Sfreamline Upwind/Pefrov-Galerkin (SUPG)- 
and Galerkin/Leasf-Squares (GLS)-sfabilized finife elemenf mefhod in combination with space-time meshes 
to discretize the full monolifhic sysfem of consfifufive equafion (3) and Navier-Sfokes equafions (2). The 
SUPG method, which has originally been proposed in [19], will serve as the stabilization method of fhe 
advection-dominafed consfifufive equafion, whereas a modified adjoint GLS will be used to stabilize 
the momentum equation [20, 21, 22]. The choice of a space-fime method is mainly motivated by future 
applications to deforming-domain problems. 

Assuming a slicing of our space-fime domain Q into N slices Q„, each spanning the computational 
domain from t„ to f„+i, and furthermore a triangulation of Q„ by fhe elements collected in we introduce 
the function space 


G C°(Q„) I VQjj e Th.m'^ ° I’qs ^ IP2 ® Pij . 

Here, the Lagrange elements P 2 and Pi are employed in space and time, respectively, with Tge denoting 
the isoparametric geometrical mapping from fhe reference elemenf onto Q[[. For all applications wifhin 
fhis paper the complete space-time domain simplifies to Q = [0, T] X Q and fhe corresponding slices fo 

Qn — [fn/ffi+l) ^ 

Furfhermore, defining fhe spatial boundary of fhe space-fime slab as ^ where Qj 

designafes fhe spatial exfenf of fhe compufafional domain af a given insfanf of time t, we use fhe following 
frial and fesf spaces 

Sh,n = {{U, p, ‘P) e X Vu,n X \ h|p„,„ = g„, 'P|p„^ = g^] (7) 

'Vh,n = {(v, q, O) e {Vh,nf X Vu,n X \ v|p„,„ = 0, Cl.|p„^ = O} , (8) 

with H and Pn,'v being the subsets of P„ on which and are prescribed as Dirichlet boundary conditions. 
The full frial space, spanning fhe whole space-fime domain, is hence chosen as 

Sh={{u,p,'¥) e L2(Q,]R‘*+i+‘b'^+i)/2) | e . 
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Using these definitions, the discretized weak problem can be formulated as follows: Given the initial 
conditions («'')“ = uq and = ‘Pq, zve seek z’' = ) e Sh such that on each time slab Q„ with 

0 < n < N ~ 1 and for every = (v^, f', e 'Vh,n the following equation is fulfilled: 

0 = an{w\z^):= J v''• p (<?,«''+ («'“• V)h'') + J ^e(v'’) : (e'^" - l) 

+ r 2Mv'’) :£(«") - r (V-v^)p>’+ r (v\Up 

dQ„ Jq„ JQn 

+ Yj f , (pCm" • V)v'' + Vq^ + psAv" - • o'*) 

■ (p(d,M'* + (m'' • V)m'') + Vp'* - psAm'' - y V • (e'^" - l)) (9) 

+ ^ p" (V •«'') + ^ (‘I*"): : ^ (('P''): - {'V^n) 

+ f^ ^(0" +w(h'*-V)0") 

: (dfT'' + (m'* • V)'P'' + [T^O(M'')] - F('P^ £(«'')) + jP{e'^") e“'*’") • 


The inner product O : T is as usual defined as fr(0^‘P). This weak form also incorporates concepts which are 
typical for space-fime GLS realizafion, e.g., fhe weak coupling between fhe space-fime slabs mofivafed by 
Disconfinuous Galerkin mefhods. Here, (h'')^ is used as fhe shorf form for lim^^o u^{tn ± <£/ •) and Q„ = Ot„. 

Furfhermore, for all subsequenf calculafions wifhin fhis paper fhe sfabilizafion paramefers were chosen 
as 


h^ h At 


'Zcons 


=min 



+ A“i 



where At is fhe fime-sfep size, h fhe elemenf diamefer, p = ps + pp fhe full viscosify, and |u| fhe absolufe 
value of the velocity evaluated at the element center. In cases where stationary simulations were performed, 
fhe corresponding parfs of fhe discrefized weak form, namely, fhe explicif fime-derivafives as well as fhe 
disconfinuous coupling across space-fime slabs, were neglecfed, which also applies fo fhe At/2 parf of fhe 
sfabilizafion consfanfs. Similarly, in fhe creeping flow limif {Re = 0) fhe advecfive derivafive of fhe velocify 
(m'' • V)m'' was omiffed from fhe momenfum equafion in conjuncfion wifh dropping h/{2\u\) from Zmom- 


3.2. Linearization and evaluation 

In a lasf sfep, fhe discrefized weak form (9) has fo be linearized in order fo make if amenable fo linear 
solvers. As already menfioned in fhe infroducfion, fhe used linearizafion mefhod in fhis work is fhe 
Newfon-Raphson algorifhm, which promises quadrafic convergence af fhe addifional cosf of providing a 
variafional direcfional derivative of fhe weak form. More specifically, denofing fhe directional derivative by 


Dan{w\-)l„ 6z|;, = — 


4=0 


a„{w^,zli + 


we iteratively solve 


Da„{w^, •) „ 6zl i = - a„{w^, z ^„,) Vw'* e ^h,n 
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for dz’lj e 'Vh.n- The updated solution can then be computed as = z*, + The iteration is 
terminated as usual when the Euclidean norm of the residual ||r ||2 := \\a„{-, z^ ;)||2 becomes smaller than a 
given threshold. 

When the Newton-Raphson algorithm is employed in the context of the newly derived constitutive 
equation (3), the immediate numerical implementation may lead to difficulties: Due to their invariance on 
the exact contour of E, the evaluation of Cauchy-type integrals is prone to cancellation. This applies to the 
double integral as well as to the exponential mapping. The difficulty can be alleviated in the numerical 
setting by evaluating the integral directly or indirectly (e.g., through a quadrature rule) only at specific 
instants of space and time. This condenses our Banach algebra to the usual matrix algebra, which essentially 
implies that the spectrum of contains at most up to d distinct discrete points, i.e., the eigenvalues of 

Using the same techniques as are applied to identify the usual spectral decomposition method of 
interpreting matrix functions with the Cauchy-type definition of matrix functions (4), we will be able to 
reformulate the integral in the framework of eigenvalues and eigenvectors. 

For this, we introduce a set of d eigenvalues A; and d eigenvectors e; of ^^(t, x), which are associated to a 
projection operator Pi = etej that projects onto the one-dimensional subspaces spanned by the corresponding 
eigenvector. Using this notation, linear algebra states that 


1 

z-T"- 


i=l 


( 10 ) 


has to hold, where for the sake of brevity the function arguments (t, x) have been dropped. Applying this 
equation to the double integral simplifies it to 




■dzdz ', 


which then together with Cauchy's integral formula (or the residue theorem) yields 

d 

F('P^ £(«")) = /(A, - Ai)P,e{u^)Pj . (11) 

w=i 


It should be noted that this is the form of the e(H)-term in the constitutive equation as it has been considered 
in [14]. 

Of course, there is no doubt that, with the numerically well-studied QR-algorithms in mind, this form is 
much more suitable for numerical evaluation. Nonetheless, it falls short in many applications when it comes 
to study perturbations of *F^, as it is the case for the variational derivative needed in the Newton-Raphson 
algorithm. Existing implementations, as, e.g., in [13,14], were limited to the two-dimensional case, since for 
a 2 X 2 matrix it is still feasible to derive an algebraic closed expression for eigenvalues and eigenvectors in 
dependence of Another approach would be general perturbation theory [23], which directly applies to 
eigenvalues A; and their projection operators P,, but this theory is prone to singularities in the vicinity of 
degenerate eigenvalues. 

The solution we will pursue here is similar to the general perturbation method in means of using complex 
calculus. As such we perform the perturbation first in the framework of the double integral and then switch 
to the eigenvalue representation. E.g., considering the variational derivative of the double integral F{'¥, e{u)) 
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with respect to '¥^ in the direction 6'¥^, one obtains by similar means as in Eq. (6) 




s=o 




' _ >ph 


_ tph 


(2m) 


1 r' r' 1 1 


-6T 


2' —'pIt 
1 


vp/i 2' — 'P^ 


dz dz'. 


Inserting Eq. (10) and applying the residue theorem then yields 




4=0 


F(t'' + 5 6'P^ £(«'')) = 

!,/,*:=1 


/(A,--A,)-/(Ay-A,) 


A, - A; 


(i’,dT'‘Pye(H'')i’, + Pte(u'^)Pj6'P'‘Pi) , (12) 


where in accordance with the residue theorem the difference quotient has to be replaced by /'(A, - Ajt) if A, 
and Ay coincide. It is clear that this formula can be evaluated along the same lines as the evaluation of the 
double integral itself (11). 

Similar considerations also yield the different derivatives of the exponential mapping as involved in the 
discretized weak form (9) 




4=0 


exp(T'' + <S6'P'‘) 




— \ A; Ay Ay A; 


Pid'P'^P 


Y-i , ,, , ,,sinh((A/- Ay)/2) 
yE/"”-"” (A.-A ,)/2 


V expCP*) = 2^ 

!,J,k=l 


(A; - A,)/2 


Here, the vectors c, denote the Cartesian basis vectors. Additionally, due to the GLS stabilization, the 
variational derivative has to be considered for V • exp('P*). The analysis yields 




4=0 


■exp(T'* + ^6T''') = ^ 


gh/2gA,/2 


sinh((A; - A,)/2) ^ ^ 


i,j,k=l 


(A, - A,)/2 


-Pidid'V^Pke, 


E 


e'kj 


p^k 


\ (A, - A;)(A,- - A,) (Ay - A,)(A; - A,) (A, - A,)(A, - \j) 

d 

■ Yj [Prb'V^Pjdi^’^Pkei + Pkdi'V^Pjb^^Pre] . 


i=i 


Eor the numerical implementation, we will have to further rewrite the factor in the second sum, as it is in 
this form not appropriate for evaluation in the proximity of degenerate eigenvalues. Introducing auxiliary 
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variables x - (A,- - y = {Ai- Ajt)/3, and z = (Ay - Ajt)/3, it can be reformulated as 


e'^i 




{A, - Aj){A^ - A,) {Aj - AdiAj - A,) (A, - A,)(A, - A;) 

1 


_^A,/3^Ay/3^At/3 


e’^-lev-l e 

-+ 


/ V'-k 'H/V'-k 


X y 


-X 


-y -z 


( 13 ) 


y-x 


/G^-1 e-y-1^ 

1 

(e^-1 e-^-1) 

1 

Uy-i e^-lV 

1 

1 

1 

x + z 

^ X -z ^ 

1 

y-z 

1 y z /. 


Thus, the evaluation is once again reduced to difference quofienfs. The laffer, which already appeared in 
Eq. (12), can be easily approximafed by a Taylor series in fhe vicinify of vanishing denominafors 


g(^) - g(y) 

x-y 



- yf 

24 




Here g{x) can sfand eifher for f{x) as in fhe case of Eq. (12) or for (e^ - l)/x as in Eq. (13). In fhe former case 
fhe Taylor series is used for \x — y\ < 10“^ and in fhe laffer case for \x - y\ < 10“^ in order fo accounf for fhe 
use of floafrng poinf arifhmefic wifh double precision. In addifion, fhe derivafives of g{x) also have fo be 
approximafed for small values of x, which for f{x) is performed beneafh |x| < 10“^ and for (e^ - l)/x below 
|x| < 10“^. The numbers were obfarned by comparing results of a naive double precision implemenfafion 
with the results computed in a much higher precision [24], and afterwards choosing the Taylor polynomials 
such that the relative error should not exceed ~ 10“^*^. A higher precision is also quite unlikely to be needed 
as both terms enter the linear equation system only on the left-hand side, such that they only influence fhe 
convergence, buf nof fhe accuracy of fhe solufion. As we will lafer see, fhe convergence is influenced even 
more by fhe inexacf solufion of fhe equafion sysfem fhrough iferafive linear solvers. 

It should be noted that such a special treatment for the other functions involved is in general unnecessary. 
The hyperbolic functions tanh(x/2) and sinh(x/2) can be readily used if j/2 + 0, assuming their implementa¬ 
tion is correctly rounded close to 0 and the rounding mode is set to round to nearest [25]. Eor {e^ - l)/j, the 
accuracy near x ^ 0 can be greatly improved by using a trick [26] and evaluating 


e^-1 

X 


y-1 

logy 


with y = , 


which is substituted by 1 in the case of y = 1. 

All other terms on the left-hand side of fhe linear equafion system arising from (9) can be derived as 
usual. The derivatives originating from the stabilization terms, in particular the derivatives of Tcons{u^ ■ V)0^' 
with respect to m*, are typically omitted, as they decrease the robustness of fhe Newfon-Raphson algorithm. 
Nonetheless, they are an essential ingredient for quadrafic convergence, which will also be discussed in fhe 
nexf section. 

In our implemenfafion we use an ILUT-precondifioned EGMRES implemenfafion fo solve fhe resulfing 
linear sysfems (cf. [27,28]). 


4. Benchmarks 

In fhis section we will use fhe newly derived mefhod to study two benchmarks: the sedimenting sphere 
benchmark and a variation thereof wifh fhe sphere replaced by a fri-axial ellipsoid. 

4.2. Sedimenting sphere 

The sedimenting sphere in a tube benchmark is, in addition to the drag on confined cylinder benchmark, 
one of the classic benchmarks that has been used in the past to measure the performance of numerical codes 
and differenf constitutive models. It has been intensively studied experimentally, as well as numerically. 
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Figure 2: Illustration of the geometry and prescribed boundary conditions for the simulation of a uniform flow past a static sphere of 
radius R. 



Figure 3: Cut through the xy-plane of Mesh Ml. 


where numerous results for the upper-convected Maxwell model were obtained. For a thorough review of 
the two aspects we refer to [29, 3] and the references therein. Our analysis will be mostly centered around 
the Oldroyd-B model, which was already analyzed in [30, 31, 32, 33[, as well as the Giesekus model. In 
contrast to the just-mentioned literature, we will not exploit the rotational symmetry in order to perform an 
in essence two-dimensional simulation of the three-dimensional problem, but will solve the problem in three 
dimensions. The latter, although computationally more expensive, is of course more flexible and preferred 
with the view on future applications. As is commonly done, we will furthermore restrict ourselves to the 
simulation of fhe fully-developed flow condifion, where fhe sphere is sedimenfing at constant speed, such 
that through a shift into the reference frame of fhe sphere, we can reformulafe fhe problem as a sfafionary 
problem of a sphere at rest within a flow with uniform velocity u. Moreover, the gravitational force is 
neglecfed; wifh fhe excepfion of a missing buoyancy ferm in fhe pressure p, fhis will nof lead to any change 
of fhe flow field. 

The geomefry, as illustrafed in Fig. 2, feafures a sphere of radius R. The sphere is locafed in fhe cenfer of 
a tube with radius 2R and is exposed to a uniform stream u = (Q, 0,0)^. Based on the geometry, the flow 
condifions, and fhe relaxafion time A, we define fhe Weissenberg number as 


It should be mentioned that the choice of the flow in x-direction is solely for the purpose of a better illustration. 
The boundary conditions, as shown in Fig. 2, are a no-slip condition on the sphere, a uniform stream of 
sfress-free polymers ('P = 0) af fhe inflow, and vanishing velocifies perpendicular fo fhe symmefry axis on 
the outflow. In accordance with the literature, only the creeping flow limif {Re — 0) is considered and the 
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Ml 

M2 

M3 

Number of elements on the sphere 

676 

2602 

9432 

Total number of nodes 

42788 

217789 

970454 

Total number of elements 

29791 

157757 

714417 

Krylov-space dimension 

150 

300 / 350 

400 

ILUT maximal fill-in nimr 

120 / 200 

120 / 200 / 250 

120 / 200 / 250 

ILUT threshold 

10-4 

10-'* 

10-4 

Number of cores 

32 

256 / 512 

2048 


Table 1: Mesh properties and solver parameters for the sedimenting sphere benchmark. 


Wi 



K 





Ml 

M2 

M3 

[30] 

[31] 

[32] 

[33] 

0.1 

5.90022 

5.90472 

5.90576 





0.2 

5.80240 

5.80646 

5.80763 





0.3 

5.68858 

5.69227 

5.69356 

5.69368 

5.6963 



0.4 

5.58068 

5.58390 

5.58527 





0.5 

5.48692 

5.48953 

5.49093 



5.4852 


0.6 

5.40899 

5.41086 

5.41227 

5.41225 

5.4117 

5.4009 


0.7 

5.34592 

5.34700 

5.34838 



5.3411 


0.8 

5.29582 

5.29616 

5.29747 



5.2945 


0.9 

5.25660 

5.25639 

5.25761 

5.25717 


5.2518 


1.0 

5.22628 

5.22586 

5.22700 



5.2240 


1.1 

5.20312 

5.20292 

5.20402 



5.2029 


1.2 

5.18568 

5.18619 

5.18733 

5.18648 


5.1842 

5.1877 

1.3 

5.17278 

5.17449 

5.17581 




5.1763 

1.4 

5.16354 

5.16689 

5.16851 





1.5 

5.15723 

5.16261 


5.15293 





Table 2: Results for the correction factor K of the drag on the sphere when using the Oldroyd-B model. 


viscosity ratio is, in all conducted simulations, chosen as jS = = 0.5. 

As already indicated in the previous section, a tetrahedral P 2 mesh was used to discretize the domain. A 
cut through the coarsest of the used meshes can be seen in Fig. 3. All meshes feature a 0.9K-thick boundary 
layer around the sphere in order to properly resolve steep gradients. Further mesh properties as well as 
the linear solver parameters can be taken from Tab. 1. Moreover, during fhe calculations fhe Weissenberg 
number was consecutively increased in such a way fhaf fhe lasf resulf always served as an initial guess for 
the following run. Run times for a single simulation range approximately from 400 s to 700 s wall-clock time 
for fhe Meshes Ml and M2 on fhe Infel-based RWTH clusfer. For fhe finesf mesh, fhe IBM-based Juqueen 
compufer was used, resulting in run times of 1300 - 2200 s. 

4.1.1. Oldroyd-B model 

The use of the described benchmark in conjunction with the Oldroyd-B model has been covered 
extensively in literature. One widely recognized performance quantity is the wall correction factor K, which 
is given as the ratio of fhe drag force on the sphere to the Stokesian drag of a sphere in an unbounded 
Newfonian fluid 


K = 


6nyRu 


f el -pi -(- 2ps£{u) -h y - 1) 

^ r Sphere 


(14) 


Here, n denofes fhe unif normal field on fhe sphere, as usual. 

The resulfs of fhe simulations as presenfed in Tab. 2 mafch fhe resulfs in liferafure quife well: Generally 
convergence can — independent of mesh size — be claimed up to a Weissenberg number of Wi = 1.4. Above 
fhis fhreshold, fhe conditioning of fhe linearized sysfem regresses. This can be mifigafed only fo a cerfain 
exfenf by an increased number of GMRES iterafions and an increased ILUT fill-in, but otherwise usually 
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Figure 4: '¥xx and o-xx plotted along the centerline in the wake of the sphere. 


leads to a failure of fhe simulafion. 

If should be nofed fhaf fhe drag mighf nof necessarily be fhe besf benchmark quanfity fo measure fhe 
performance of numerical discrehzafions, which may yield perfecf drag resulfs while not being able to 
properly predict other important flow characferisfics. One of fhese characterisfics is fhe exfensional flow 
in the wake of fhe sphere, where fhe polymers are strefched along fhe flow direction. This is of special 
importance for the Oldroyd-B model, which, as it corresponds to the microscopic Hookean-dumbbell model, 
has the property that the solution blows up in a purely extensional flow if the extensional rate exceeds a 
critical point — simply put, the dumbbells become infinitely long. Although there has not been a conclusive 
proof in liferafure yef, if is believed fhaf a similar mechanism is also responsible for fhe limifafion in fhe 
Weissenberg number for the feasible simulations in the falling sphere benchmark. To highlight this similarity, 
notice that symmetry dictates T and Vh to be diagonal along the centerline, such that the constitutive 
equation of 'Vxx in Eq. (3) reduces fo 


Uxdx'i'xx - 2dxUx = (l - e . (15) 

Considering that at any extremal point x* of '¥xx the derivative has to vanish, rearranging this equation yields 

'i'xxix*) = - log (1 - 2A dxUxix*)) . 

Thus, with dxUx{x*) approaching 1/(2A) the component 'Vxx blows up. Of course, nofhing parficular 
on fhe behavior of dxUx can be inferred within the framework of one-dimensional analysis due fo fhe 
incompressibilify consfrainf. 

The 'Vxx acfually predicfed by fhe simulafion can be seen in Fig. 4. One of fhe poinfs fhat becomes 
directly apparent is that in these simulations mesh convergence can only be claimed up to Weissenberg 
numbers Wi = 1.0 - 1.2. Above these values, it seems that despite a boundary-layer-resolvmg mesh, the 
fluid characferisfics in fhaf region sfill cannof be accurafely described. This effecf becomes even more 
pronounced looking af (Txx, which modulo numerical noise is given by (Txx = exp 'Vxx arid is also depicfed in 
Fig. 4. There, of course, a slight deviation of an already large 'Vxx is furfher amplified by fhe exponenfial 
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Wi 





K 






a = 0.001 



a = 0.01 



a = 0.1 


Ml 

M2 

M3 

Ml 

M2 

M3 

Ml 

M2 

M3 

0.1 

5.89918 

5.90369 

5.90473 

5.88997 

5.89464 

5.89573 

5.81454 

5.82032 

5.82166 

0.2 

5.79863 

5.80274 

5.80393 

5.76691 

5.77147 

5.77275 

5.56351 

5.57002 

5.57160 

0.3 

5.68098 

5.68479 

5.68610 

5.62095 

5.62552 

5.62694 

5.31523 

5.32188 

5.32349 

0.4 

5.56845 

5.57185 

5.57324 

5.47847 

5.48300 

5.48451 

5.09969 

5.10625 

5.10785 

0.5 

5.46934 

5.47222 

5.47366 

5.34928 

5.35374 

5.35531 

4.91688 

4.92331 

4.92489 

0.6 

5.38538 

5.38763 

5.38910 

5.23534 

5.23965 

5.24127 

4.76150 

4.76781 

4.76938 

0.7 

5.31561 

5.31714 

5.31861 

5.13538 

5.13952 

5.14118 

4.62842 

4.63461 

4.63616 

0.8 

5.25811 

5.25891 

5.26037 

5.04716 

5.05109 

5.05280 

4.51345 

4.51955 

4.52109 

0.9 

5.21075 

5.21091 

5.21235 

4.96839 

4.97211 

4.97387 

4.41335 

4.41936 

4.42090 

1.0 

5.17150 

5.17122 

5.17264 

4.89714 

4.90066 

4.90248 

4.32556 

4.33150 

4.33303 

1.1 

5.13857 

5.13811 

5.13955 

4.83192 

4.83528 

4.83716 

4.24805 

4.25394 

4.25545 

1.2 

5.11050 

5.11012 

5.11165 

4.77169 

4.77491 

4.77684 

4.17920 

4.18503 

4.18653 

1.3 

5.08611 

5.08608 

5.08774 

4.71569 

4.71879 

4.72077 

4.11769 

4.12347 

4.12496 

1.4 

5.06448 

5.06501 

5.06688 

4.66336 

4.66637 

4.66840 

4.06245 

4.06818 

4.06966 

1.5 

5.04489 

5.04617 

5.04829 

4.61432 

4.61725 

4.61931 

4.01260 

4.01828 

4.01975 

1.6 


5.02897 

5.03139 

4.56825 

4.57111 

4.57319 

3.96740 

3.97303 

3.97448 

1.8 




4.48402 

4.48676 

4.48886 

3.88863 

3.89413 

3.89557 

2.0 




4.40900 

4.41166 

4.41375 

3.82234 

3.82771 

3.82914 

2.2 




4.34189 

4.34447 

4.34653 

3.76580 

3.77103 

3.77245 

2.4 




4.28156 

4.28409 

4.28612 

3.71703 

3.72211 

3.72352 

2.6 




4.22709 

4.22960 

4.23159 

3.67452 

3.67943 

3.68085 

3.0 




4.13269 

4.13526 

4.13716 

3.60400 

3.60857 

3.61001 

3.5 




4.03595 

4.03881 

4.04057 

3.53552 

3.53969 

3.54117 

4.0 




3.95673 

3.96008 

3.96170 

3.48206 

3.48589 

3.48740 

4.5 




3.89055 

3.89454 

3.89602 

3.43910 

3.44264 

3.44418 

5.0 




3.83433 

3.83905 

3.84040 

3.40376 

3.40707 

3.40864 

5.5 




3.78591 

3.79140 

3.79265 

3.37414 

3.37728 

3.37886 

6.5 




3.70655 

3.71357 

3.71469 

3.32716 

3.33010 

3.33166 

7.5 




3.64403 

3.65246 

3.65355 

3.29144 

3.29432 

3.29584 

8.0 




3.61743 

3.62651 

3.62762 

3.27658 

3.27947 

3.28096 

8.5 




3.59334 

3.60301 

3.60415 

3.26328 

3.26619 

3.26765 

9.0 




3.57140 

3.58160 

3.58280 

3.25130 

3.25424 

3.25567 

9.5 




3.55134 

3.56202 

3.56328 

3.24044 

3.24344 

3.24482 

10.0 




3.53290 

3.54402 

3.54535 

3.23055 

3.23360 

3.23495 

11.0 




3.50015 


3.51352 

3.21318 

3.21637 

3.21764 

12.0 




3.47192 


3.48607 

3.19839 

3.20173 

3.20293 

13.0 




3.44729 


3.46211 

3.18564 

3.18914 

3.19027 

14.0 




3.42558 


3.44098 

3.17452 

3.17816 

3.17923 

15.0 




3.40628 


3.42218 

3.16472 

3.16851 

3.16952 


Table 3: Results for the correction factor K of the drag on the sphere when using the Giesekus model. 


function. Furthermore, it should also be noted that the point where attains its maximum seems to reach 
its maximal x-value at Wi - 0.8. For higher Weissenberg numbers the maximum is then shifted again in the 
direction of the sphere. 

Considering these peculiarities and the general non-linearity of the governing equations, it is even more 
remarkable that on the downward slope of in Fig. 4 the field shows a linear behavior. Performing a least 
squares fit of linear curves to the simulated data in the region x = lOP - 12R on Mesh M3 yields slopes of 
m = -0.6259 for Wi = 1.4, m = —0.7037R~^ for Wi = 1.2, and m = -0.7708 for Wi = 1.0. The resulting 
linear curves can be examined in Fig. 4. It is yet unclear which mechanism leads to this linear behavior. 

4.1.2. Giesekus model 

A similar analysis as for the Oldroyd-B model is also conducted for the Giesekus model that extends 
the Oldroyd-B model by an additional term. In fact, the Oldroyd-B model is a special case of the Giesekus 
model for a vanishing mobility parameter a = 0. 
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Figure 5: Wall correction factor K plotted for different values of the mobility a, computed on the Mesh M3. 


As before, the quantity studied first is the drag correction factor K for several Weissenberg numbers 
Wi and varying mobility a. The results are collected in Tab. 3 and depicted in Fig. 5. It is notable that, for 
the smallest a = 0.001, the model exhibits similar numerical behavior as the Oldroyd-B model, namely 
convergent results only up to Wi — 1.6. This is to be expected, higher Weissenberg numbers may be 
achievable with finer meshes in contrast to the Oldroyd-B model. Increasing a shows that the drag on 
the sphere decreases in general, which is attributable to the shear-thinning properties of a Giesekus fluid. 
Moreover, for all performed numerical calculations the drag is monotonically decreasing with increasing 
Weissenberg number Wi, and there are indications that K reaches a plateau for sufficiently high Wi. 

Looking at the extensional flow characteristics of the Giesekus model in the wake of the sphere, it already 


(a) Wi = 2.0 'Vxx 


m 


B O.85 

-0.40 

- 1.66 

5.42 

3.24 

1.05 

-1.13 

-3.31 

Figure 6: Cut through the ry-plane of Mesh M3, illustrating Wxx for different Weissenberg numbers and a = 0.1. 
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Figure 7: *Pyy and Ux plotted along the centerline in the wake of the sphere for a = 0.1. 


becomes apparent from the equivalent of Eq. (15) that the model is better-behaved: 

Uxdx'Vxx - '2-dxUx = -^ (l - 2a - (1 - a)e~'*° + ae'^“) . 

Here, the additional a exp(4'j;j:) term can potentially compensate an increase of dxUx exceeding 1 /(2A), thus 
limiting the increase of '¥xx- The resulting computations of '¥xx for two different Weissenberg numbers 
Wi = 2.0 and Wi = 15.0 are shown in Fig. 6. The results reflect clearly that with increasing Weissenberg 
number, the polymers need more time to relax to their stress-free state, which means that they are transported 
further downstream before they reach this state. As such, the 'Vxx contours also extent further downstream 
for higher Weissenberg numbers than for lower ones. As a consequence the demands on the used geometry 
and meshes increase: They need to sustain a high refinement level over a larger region in the wake of the 
sphere. 

This effect becomes even more noticeable when considering the other degrees of freedom in our simulation. 
In Fig. 7, *Pyy has been plotted along the centerline for different Weissenberg numbers. The first point to 
notice is that mesh convergence can be reached within the boundary-layer-adjusted mesh around the sphere, 
but as soon as the mesh resolution decreases, the accuracy in the to-be-predicted degree of freedom *Pyy is 
lost. The impact becomes more severe the higher the Weissenberg number is. In addition, by inspecting 
Fig. 7, it seems that for Wi = 15.0, Tyy exhibits a small kink around x = 10.5 R on Mesh M3, which may be 
attributable to a still insufficient refinement level of the mesh in that particular region. 

On the other hand, the fact that Tyy is negative also means that errors therein are exponentially damped 
in (Tyy = exp *Pyy. Smce the latter is what essentially contributes to the momentum equation, it is not much 
of a surprise that the velocity component depicted in Fig. 7 is still smooth for all used meshes. Furthermore, 
velocity overshoots exceeding u are clearly visible, in contrast to Oldroyd-B simulations. Moreover, the 
downstream relaxation is once more delayed with increasing Weissenberg number. 

4.1.3. Performance of the Nezvton-Raphson algorithm 

It is important to note that the Newton-Raphson method indeed delivers quadratic convergence up to 
the point where the errors emerging from the inexact linear solution or the limited floating-point accuracy 
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(a) m = 0.3 and a = 0.001 


(b) Mesh Ml and a = 0.001 




Newton-Raphson iteration Newton-Raphson iteration 

Figure 8: Convergence behavior of the Newton-Raphson algorithm for different settings. 


become dominant. Fig. 8 depicts a comparative study of the residual after each Newton-Raphson iteration 
across different meshes, as well as an analysis of the convergence behavior for different Weissenberg numbers. 
In both cases, the residual has been evaluated in the Euclidean norm and scaled by the square root of fhe 
total number of degrees of freedom in order to make the results comparable across different mesh sizes. 
As one sees in Fig. 8a, the convergence history can be roughly split into three phases: In the first step, 
the improvement is rather moderate; at most linear convergence was obtained. In the second phase, e.g., 
for the Mesh M3, one observes a relative improvement of the residual by a factor of 18.34 in the second 
step and an improvement by a factor of 1097 > (18.34)^ in fhe fhird sfep; fhus slighfly exceeding quadratic 
convergence. If becomes apparenf from Fig. 8a fhaf fhis convergence is mesh-independenf. The fhird phase 
is fhen domrnafed by errors infroduced by fhe inexacf solution of fhe linear equation sysfems. This can be 
deduced from the results in Fig. 8a, where the calculations on Mesh M2 were performed using two different 
ILUT fill-in settings. A reducfion of fhe ILUT fill-in, and a fherewifh increased error in fhe solution of fhe 
linear systems, directly leads to a deterioration of the quadratic convergence to at most asymptotically linear 
convergence [34]. Furthermore, in the last steps the convergence is limited by the fact that a residual far 
beneafh 10“^® is in general nof attainable due fo fhe floafing-point arifhmefic used. 

In Fig. 8b, one notices that, using the result obtained for the previously calculated Weissenberg number as 
an initial guess for the subsequent calculation, the convergence progression is similar across the consecutive 
runs. Only the starting point Wf = 0.1 does not fully fit into this picture, which on the one hand has to 
be attributed to the circumstance that for this case the initial guess was set to zero in the interior of the 
computational domain, and on the other hand is a consequence of the derivative of TconsW ^' V)®^' in the 
discretized weak form being neglected (cf. Section 3.2). The latter is a remedy for the fact that without these 
additional terms, the iterative scheme seems to be more robust with regard to the choice of fhe initial guess. 

4.2. Sedimenting ellipsoid 

In order to demonstrate the applicability of the proposed method to a truly three-dimensional problem, a 
case similar to the sedimenting sphere benchmark is considered, but with the sphere replaced by a tri-axial 
ellipsoid. The latter was chosen with the semi-principal axes aligned to the coordinate axes. The length 
of the axes in x,y, and z direction were set to a = 1.25 R,b = 1.0 R, and c = 0.8 R, respectively, such that in 
accordance with the findings in [35], fhe semi-major axis coincides wifh fhe main flow direction for Re = 0. 
The fube radius was kepf as 2R. 
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Figure 9: Cut through the xy-plane of Mesh M4, as used in the calculations with an ellipsoid. 


M4 

Number of elements on the ellipsoid 

2706 

Total number of nodes 

218663 

Total number of elements 

158353 

Krylov-space dimension 

300 

ILUT maximal fill-in h/lut 

200 

ILUT threshold 

10“^ 

Number of cores 

256 


Table 4: Mesh and solver attributes used for the sedimenting ellipsoid benchmark. 


The mesh, as depicted in Fig. 9, was chosen similar to the Mesh M2 in the sedimenting sphere benchmark, 
which already provided a good trade-off between computational cost and accuracy in the drag computation. 
Therefore, the GMRES/ILUT parameters were also chosen accordingly, as can be seen in Tab.4. 

Our main objective of the investigation was the drag correction factor K, where the latter has been 
defined for the sake of simplicity as in the case of the falling sphere, cf. Eq. (14). Nonetheless, the Stokesian 
drag formula can be generalized to ellipsoids in principle [36]. The results in Tab. 5 confirm the general 
trend of the simulations with the spherical geometry: The drag decreases monotonically with increasing 
Weissenberg number. It can also be stated that the general drag level is below the drag levels obtained 
in the simulations with a sphere as obstacle, which may be attributed to the reduced cross section. With 
increasing a, higher Weissenberg numbers can be attained, and the effect of reduced drag due to increased 
shear-thinning becomes visible. 

5. Conclusion and discussion 

The main objective of this paper was to derive a log-conformation formulation that on the one hand 
inherits the stability properties of the originally proposed log-conformation formulation [8], but on the other 
hand also paves the way for an application of Newton's method in numerical simulations. Furthermore, 
we especially sought a description that could be applied in three dimensions with the same ease as the 
previously published two-dimensional approaches [13,14]. 

To demonstrate the numerical benefit of this approach, we implemented a proof-of-concept three- 
dimensional finite element solver and subsequently tested it by means of the sedimenting sphere and 
ellipsoid benchmarks. The simulations exhibited the best-possible convergence properties of quadratic- 
convergence. 

Since the new constitutive equations are just a rewording of the original log-conformation equations, 
the proposed formulation cannot further improve the stability. As such, we were not able to obtain results 
beyond a Weissenberg number of Wi = 1.4 for a sphere sedimenting through an Oldroyd-B fluid. Since 
switching to the Giesekus model removed this limitation, the characteristic behavior of the Oldroyd-B fluid 
in extensional flow regimes might be the underlying reason for this restriction. 

In addition to the just-mentioned advantages for the numerical application, our formulation is intrinsically 
defined in an imdiscretized setting, which may reveal new perspectives on the analytical properties of the 
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Wi 


K 



a = 0 

a = 0.001 

a = 0.01 

a = 0.1 

0.1 

4.90847 

4.90782 

4.90211 

4.85331 

0.2 

4.85959 

4.85715 

4.83621 

4.68884 

0.3 

4.79819 

4.79305 

4.75092 

4.50990 

0.4 

4.73648 

4.72791 

4.66130 

4.34532 

0.5 

4.68082 

4.66810 

4.57513 

4.20056 

0.6 

4.63365 

4.61608 

4.49540 

4.07438 

0.7 

4.59537 

4.57217 

4.42252 

3.96425 

0.8 

4.56535 

4.53566 

4.35581 

3.86772 

0.9 

4.54252 

4.50538 

4.29432 

3.78269 

1.0 

4.52568 

4.48007 

4.23718 

3.70739 

1.1 

4.51370 

4.45854 

4.18372 

3.64038 

1.2 

4.50554 

4.43976 

4.13345 

3.58045 

1.3 

4.50034 

4.42286 

4.08604 

3.52658 

1.4 

4.49738 

4.40717 

4.04124 

3.47795 

1.5 


4.39215 

3.99883 

3.43385 

1.6 


4.37747 

3.95867 

3.39369 

1.8 



3.88450 

3.32332 

2.0 



3.81772 

3.26370 


10.0 

3.01520 

2.71617 

11.0 

2.98506 

2.69985 

12.0 

2.95902 

2.68597 

13.0 

2.93626 

2.67400 

14.0 

2.91617 

2.66356 

15.0 

2.89827 

2.65437 


Table 5: Results for the correction factor K of the drag on the ellipsoid. 


used constitutive models in the future. In particular, the seamless incorporation of fhe so-called free-energy 
esfimafes in fhe log-conformafion formulafion [37] and fheir application fo fhe global-in-fime exisfence of 
solutions [38, 39] may give new insighfs. 
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